Numerical Replica Limit for the Density Correlation of the Random Dirac Fermion 
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The zero mode wave function of a massless Dirac fermion 
in the presence of a random gauge field is studied. The den- 
sity correlation function is calculated numerically and found 
to exhibit power law in the weak randomness with the disor- 
der dependent exponent. It deviates from the power law and 
the disorder dependence becomes frozen in the strong ran- 
domness. A classical statistical system is employed through 
the replica trick to interpret the results and the direct eval- 
uation of the replica limit is demonstrated numerically. The 
analytic expression of the correlation function and the free 
energy are also discussed with the replica symmetry breaking 
and the Liouville field theory. 
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Although the scaling theory of localization for two- 
dimensional disordered systems generally predicts the ab- 
sence of extended states, we have some examples of non- 
localized states in two dimensions which are marginally 
allowed to appear. Among them, systems with chiral 
symmetry such as the Anderson bond-disordered model 
|y, the random flux model or the 7r-fl.ux model 

with link disorders |||| have attracted lots of attention. 
Density of states of the models have singularities at the 
zero energy and the corresponding wave functions ex- 
hibit multifractal behavior. Actually, one generally ex- 
pects the existence of zero energy states for these models. 
Consider a Hamiltonian of the above models with chiral 
symmetry TL = X)<«> c l^ij c j + ^- c - on a bipartite lattice 
A which can be decomposed into two sublattices A^ and 
Ab- After performing a unitary transformation (redefi- 
nition of indices) , the Hamiltonian is expressed as 



w = ({<£} {4}) 



A 

AB 



V AB 
Ob 



{ca} 
{cb} 



The off-diagonal block structure of TL implements the fact 
that hopping is restricted between the interpenetrating 
sublattices A^ and A^. The zero mode wave functions 
tp = t (ipA,ipB) satisfy the Schrodinger equations 



Vab^b = 0, 



0. 



(1) 



Let Na,b is the number of sites on Aa,b- For the cases 
where Na — Nb > ( which can be realized, for example, 
by an appropriate boundary condition ), standard linear 
algebra tells us there always exist (Na — Nb) independent 
zero energy solutions for Eq. (Q) with vanishing ips- In 
this expression, the notion of chiral symmetry is explicit. 



However, one needs to solve Eq. (Q) numerically to go 
further. 
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FIG. 1. Typical probability densities |V>(:c)| 2 for g — 0.4 
on a (a) 64 x 64 lattice and (b) 1024 x 1024 lattice. Insets: 
|i/>(:r)| 2 for g = 6.4 for each system. 



Now, let us concentrate our attention on cases where 
the low- lying physics is described by the Dirac fermions. 
Remarkably, in this situation, the explicit construction 
of zero energy wave functions is possible in a two- 
dimensional continuum space Let us consider a 
Hamiltonian of the form H = a ■ p + a ■ A where <Ji =x ,y 
are the usual 2x2 Pauli matrices and A is a random 
gauge field. The Schrodinger (Dirac) equations for the 
zero modes are 



(2d + iA)iP- 
where d = {d x 



0. 
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id y )/2, A 







iA y and A 



ldy)/2, 

A x — iA y . This is a continuum analogue 
of Eq. (Q). If we adopt the Coulomb gauge to express 
the vector potential A in terms of a scalar potential $ 
as A x = d y Q, A y = —d x Q and assume that the mean 
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total flux piercing the system is zero, we can obtain the 
exact solution for any realization of disorders as ip± (x) = 
C ± e=F*( x ). Further, let us assume the probability weight 
for each realization of <I>(x) has the form P[$] oc e~ s 
where S[<&] — l/2g J d 2 x(V<I>(x)) 2 and g is the disorder 
strength or, in the field theoretic language, the coupling 
constant, which is dimensionless in two dimensions. 

From now on, we concentrate our attention only on 
i/> + . For physical interest, it is necessary to consider 
normalized wave functions in a L x L box as ^(x) = 
e -*W/VZ with Z = J^ e - 2 *( x ) where a is a lat- 
tice constant. Here, we regularized the problem on 
a N x N periodic lattice although the original prob- 
lem is formulated in a continuum space (where L = 
Na) . Correspondingly, we use the probability weight 
5[$] oc l/2g , X^<i J >(^i — ®j) 2 or > m tnc momentum 
space, S[$] cx iV 2 /.9£ m $ m $_ m (2 - Ej=i 

where $ m = iV~ 2 £\ e~ lkmX: > and the sum extends 
over the first Brillouin zone rn M = —N/2 + 1 • • • N/2 with 
kf/ n = 2mn fJ -/aN (N is even for convenience.) ||. Typical 
probability densities |?/;(x)| 2 calculated numerically are 
shown in Fig. [l], which remind us of multifractal states 
found at a localization-delocalization transition for sev- 
eral systems Q. 

In fact, the multifractal property of this wave func- 
tion has been revealed quantitatively by a close analogy 
to a generalized random energy model (l^JT^]. As the 
disorder strength g varied, the multifractal spectrum ex- 
hibits a sharp transition which is similar to the freezing 
phenomenon in spin glasses. Several other approaches 
such as the supersymmetry (SUSY) technique |13) , the 
connection to the Liouville field theory H] , the renor- 
malization group (RG) jlfj or conformal field theory filf | 
have also been taken to support the transition. 

Since the calculated probability densities (Fig. [j]) are 
so spiky, the discretization procedure above may not be 
justified. In spite of this subtlety, we concentrate on this 
well defined discretized wave function to investigate uni- 
versal properties. 

In this letter, we evaluate the density correlation func- 
tion 



<^ 2 (1)^ 2 (2)) = ^I e -»H- 



2*(x 2 



where (• • •) denotes the averaging with respect to the 
weight -P[$]. Here, the difficulties reside in the normal- 
ization factor Z in the denominator since Z itself is a 
random variable. The one of the simplest attempts to 
cope with it is the replica trick. We multiply the numer- 
ator by Z n and consider 

ty a (l)V a (2)> n = 

-2[$(xi)+*(x 2 )+*(£i)+-+*(£ n _ 2 )] 



d 2 ^ d 2 C„- 



(2) 



use this replica trick to interpret the direct numerical re- 
sults and also try to take the replica limit by evaluating 
(^ 2 (1)^ 2 (2)) numerically for several n and extrapolat- 
ing them to n — 0. In addition, we utilize the evaluation 
of (V 2 (1)V' 2 (2)) J1 together with the Liouville field theory 
to get the analytic expression of the correlation function 
for the weak disorder regime. 
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which is expected to reduce to (i/ ,2 (l)V' 2 (2)) by taking 
the replica limit n — > (analytic continuation). We 



|x, 2 |/a 

FIG. 2. The numerically calculated density correlation 
for (a) the weak disorder regime and (b) the strong disorder 
regime on a 1024 x 1024 lattice. For g < 1.6, the statistical 
error is smaller than the data point. For g > 3.2, the error 
bars are shown for every 100 points. Insets: The "exponent" 
Ajv v.s. 1/N for (a) g = 0.4 and (b) g — 6.4. It is evaluated 
for 1 < |xi2|/a < 10 for each finite system. 



First, let us present the direct numerical calculation 
of (V' 2 (1)V' 2 (2))- In this problem, the probability weight 
itself is diagonal in the momentum space (see above), 
which allows us to carry out numerical simulations with 
a very large lattice up to 2048 x 2048. Fig. || shows calcu- 
lated (ip 2 {l)ip 2 {2)) for various g on a 1024 x 1024 lattice. 
The quenched averaging is performed over ~ 10 5 differ- 
ent realization of disorders |l7|] . As is shown, the corre- 
lation function for the weak disorder exhibits power law 
behavior (<0 2 (1)^ 2 (2)) - |x 12 |" A for 1 < \x l2 \/a < A^/2 
with its exponent A dependent linearly on g, A = 2g/ir 
(Fig. H). It is consistent with the several analytic ap- 
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proaches As g increases, however, the g dependence 
of the correlation function becomes weaker and it devi- 
ates from the power law. To be more precise, there is a 
systematic deviation from the simple power law, that is, 
if we determine the "exponent" An on a finite N sys- 
tem, it seems to diverge as N increases ( see the insets 
of Fig. 2 ). It is clearly different from the behavior of 
Ajv in the weak randomness where A^v seems to con- 
verge. In fact, as is shown in Fig. [I], the wave function 
becomes peaked on few sites as g increases. However, it 
is different from the usual localized wave function which 
decays exponentially with its typical length scale char- 
acterized by the localization length. The above change 
of behavior in the correlation function is consistent with 
the transition from the weak to strong disorder found in 
the multifractal spectrum by the previous studies. 



(ip 2 (l)ip 2 (2)) n can be interpreted as the two body den- 
sity of a classical statistical system consisting of a set of 
particles (replicas) interacting each other via the poten- 
tial G{xi,Xj) fll8| . These replica estimations are shown 
and directly compared to the numerical results (Fig. ||). 
In Fig. |J, (i/) 2 (l)^ 2 (2)) n for various n (the number of 
replicas) with fixed g are obtained by calculating Eq. (j^) 
numerically on a 64 x 64 lattice. This rather small lattice 
size is due to the multiple integral in Eq. (||). Note that 
we do not have (^(l)^ 2 ^)) as is inferred from Eq. 
(^). For g = 0.4, the replica estimation seems to converge 
to the one calculated by the direct numerical simulations. 
For g = 6.4, in contrast, it hardly seems to coincide to 
the exact one in the limit n — ► 0. Moreover, it gives an 
unphysical result, i.e., a negative exponent, after taking 
the replica limit. 




FIG. 3. The "exponent" Ajv of the two-point correlation 
function with respect to the disorder strength g and its naive 
replica estimates A n =2,3,4,5 (log-log plot). The "exponent" is 
evaluated for 1 < \x 12 \/a < 10 for 64 x 64 (o), 256 x 256 (O) 
and 1024 x 1024 (•) system. The statistical error is smaller 
than the symbol size. The analytic estimation Eq. (||) for A n 
is represented by lines. 



Next, let us try to interpret the above numerical results 
by the replica trick. After replicating Z, we can perform 
the averaging (• • •) in Eq. (^) and obtain 



<V 2 (1)<A 2 (2)) „ =<*(&-ziW&-s a )>? 



(3) 



with 



(0)t 



Tr = 



Tr [Oe- H "\ /Tr [e^] 

1 fcP^ d 2 Cn 



k<l 



where G(xi,Xj) = (&(xi)&(xj)) is the Green's function 
and Q[xi,Xj) = G(xi,Xj) — G(0) ~ — <?/27rm (\xij\/a). 
Here, we multiply some trivial factors which reduce to 
unity in the limit n — ► 0. As is suggested in Eq. (0), 
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FIG. 4. The numerically calculated density correlation 
(•) and their replica estimates with the different number of 
replicas n = 2, 3, 4, 5 for (a) g = 0.4 and (b) g — 6.4 on a 
64 x 64 periodic lattice (log-log plot). The replica estimates 
are obtained by evaluating the multiple integral in Eq. (^) di- 
rectly. (^p 2 (xi)^ 2 (22)) _ Q is obtained by extrapolating from 



Insets: (^jj 2 (xi)'ip' 2 : (x2)) v.s. n at 



{^(sOV^))^,,, 

I £12 1 /a = 20 (semi- log plot) 



This breakdown is closely related to the transition in 
the replica space. There are two distinct phases for this 
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system. For small g, all configurations are equally favor- 
able. As g increases, however, the configurations where 
all replicas are close to each other come to have large 
weight. Thus, for sufficiently small g, it is enough to 
concentrate on only and in Eq. (pM. The configura- 
tion of the other particles is smeared out in the ensemble 
average and irrelevant for (^> 2 (1)^ 2 (2)) . For large g, 
on the other hand, the main contributions in the ensem- 
ble average are from the configurations where all replicas 
except (or are at x\ ( or 12, respectively). Then, 
(V' 2 (l)V' 2 (2)) n is expected to behave as 



with 



An 



J 2 g/n for small g 

\ 2g(n — I) /it for large g 



(4) 



by using the replica estimates. We expand e~^ z to ex- 
press (V' 2 (l)'0 2 (2)) as the superposition of the replica 
estimates with a different number of replicas, i.e., the 
grand canonical ensemble 

/>oo 

<^ 2 (1)^ 2 (2)) = 2 / dji 



£(-i)» e -*.oo / J2 5{tk - - x 2 ) ) 

n=2 \k^l I n 

where F n {Ji) = - \ n T re -(H n -^n) + 2 n 2 G(0) and n = e». 
Since for the weak disorder, the summand only trivially 
depends on n, we can easily sum up this suggestive ex- 
pression and obtain the same result as the replica trick. 

This method is applicable also for the free energy. We 
expand In Z = f™ & (e - " - e 



as 



Numerically calculated A n is shown in Fig. ||, which 
confirms the above estimation for small n is reasonable. 
Taking the replica limit n — ► 0, we obtain A n= o — 2g/ir 
for small g which is consistent with the results obtained 
by SUSY technique fl3f| - For the strong disorder regime, 
however, the exponent reduces to —2g/i: which is un- 
physical. 

Does it mean the replica trick is a mere trick ? One of 
the possible scenarios is that the apparent breakdown 
is because we only take the small number of replicas 
into account. So, for large n, the n dependence of A„ 
may deviate from Eq. (|4|) and give the correct answer 
in the replica limit J19| . Moreover, the replica symme- 
try breaking (RSB) solution which was proposed for the 
free energy ]15| may be applicable also for the correlation 
function. 

We also investigated other types of correlation func- 
tions such as (-0(1)^(2)) or ($(1)V> 2 (1)$(2)V> 2 (2)) , the 
latter of which is of interest because it is related to 
the second derivative of the free energy (In Z) / In (L /a) 
which shows the non-analyticity at g = 2it. Their be- 
haviors are qualitatively similar to that of (?/> 2 (1)t/> 2 (2)) 
in that, for small g, these correlation functions become 
steeper and steeper as g increases whose g dependence 
are calculable by the replica trick. For large g, how- 
ever, their g dependencies are rather weak and the naive 
replica trick fails. 

Another interesting approach is to utilize the formula 



1 

z N ~ (N-iy. 
tion function as 



J °° dfj,e ^ z fi N 1 and express the correla- 



(<A 2 (i)^(2)) 



-2*(xi)-2*0 2 ) -s L 



-[*] 



where S 



LFT 



I 



d 2 i 



-2*(i) 



It 



This 



l/2 ff (V$r + /ie- 

action resembles that of the Liouville field theory in two- 
dimensional quantum gravity. However, since it was 
pointed out that there are some subtleties about the field 
theoretic treatment of Eq. (^) , we evaluate it directly 



(hxZ) 



djl 



E( 

n=l 



■l)"e 



It is difficult to perform the summation for the strong 
disorder though we can obtain the correct answer for 
the weak disorder. However, if we simply employ the 
RSB estimate by Carpentier and Doussal fla 



1 (L/a) 



{pg/-n+2/p+p,)n 



where p = 1 for the weak disor- 



der and \J2-k / 1 g for the strong disorder, the summation 
reproduces the exact result both for the weak and strong 
disorder regime fl9[| . 
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